
# Quantile regression for IPL-12 ------------------------------------------

# Predict quantile data
taus <- c(.05,0.10, 0.25, 0.5, 0.75, 0.90,.95)
nvars <- ncol(df)
for (i in 1:length(taus)) {
  
  qreg  <- rq(ipl12s ~ yearsIn + I(yearsIn^2) + I(yearsIn^3), data = df, tau = taus[i])
  qregFit <- predict(qreg, newdata = df)
  df <- cbind(df, qregFit)
  colnames(df)[nvars + i] <- paste("per", taus[i], sep = "")
}

QRplot12 <- ggplot(df, aes(x = yearsIn, y = ipl12s)) +
  geom_jitter(aes(), size = .2, colour = 'grey') +
  geom_line(aes(x = yearsIn, y = per0.05), colour = '#8da0cb', linetype = "longdash") +
  geom_line(aes(x = yearsIn, y = per0.1), colour = '#8da0cb', linetype = "longdash") +
  geom_line(aes(x = yearsIn, y = per0.25), colour = '#8da0cb', linetype = "longdash") +
  geom_line(aes(x = yearsIn, y = per0.5), colour = '#fc8d62') +
  geom_line(aes(x = yearsIn, y = per0.75), colour = '#8da0cb', linetype = "longdash") +
  geom_line(aes(x = yearsIn, y = per0.9), colour = '#8da0cb', linetype = "longdash") +
  geom_line(aes(x = yearsIn, y = per0.95), colour = '#8da0cb', linetype = "longdash") +
  # geom_smooth(method = "loess", size = .5) +
  coord_cartesian(ylim = c(0, 1)) +
  ylab("IPL-12 Integration Index") +
  xlab("Years of Residency") +
  ggtitle("") +
  theme_minimal() +
  theme(legend.position = "none",
        axis.title = element_text(size = 8),
        axis.text = element_text(size = 7),
        plot.margin = margin(5.5, 5.5, 5.5, 22, "pt")) 

# Difference between 1st, 20th and 40th year of residence
mean(df$per0.5[df$yearsIn == 1])
mean(df$per0.5[df$yearsIn == 20])
mean(df$per0.5[df$yearsIn == 40])
mean(df$per0.5[df$yearsIn == 40]) - mean(df$per0.5[df$yearsIn == 20])
mean(df$per0.5[df$yearsIn == 20]) - mean(df$per0.5[df$yearsIn == 1])

# Quantile regression for IPL-24 ------------------------------------------

nvars <- ncol(df)
for (i in 1:length(taus)) {
  
  qreg  <- rq(ipl24s ~ yearsIn + I(yearsIn^2) + I(yearsIn^3), data = df, tau = taus[i])
  qregFit <- predict(qreg, newdata = df)
  # assign(paste("per", taus[i], sep = ""), qregFit)
  df<- cbind(df, qregFit)
  colnames(df)[nvars + i] <- paste("per24", taus[i], sep = "")
}

QRplot24 <- ggplot(df, aes(x = yearsIn, y = ipl24s)) +
  geom_jitter(aes(), size = .2, colour = 'grey') +
  geom_line(aes(x = yearsIn, y = per240.05), colour = '#8da0cb', linetype = "longdash") +
  geom_line(aes(x = yearsIn, y = per240.1), colour = '#8da0cb', linetype = "longdash") +
  geom_line(aes(x = yearsIn, y = per240.25), colour = '#8da0cb', linetype = "longdash") +
  geom_line(aes(x = yearsIn, y = per240.5), colour = '#fc8d62') +
  geom_line(aes(x = yearsIn, y = per240.75), colour = '#8da0cb', linetype = "longdash") +
  geom_line(aes(x = yearsIn, y = per240.9), colour = '#8da0cb', linetype = "longdash") +
  geom_line(aes(x = yearsIn, y = per240.95), colour = '#8da0cb', linetype = "longdash") +
  # geom_smooth(method = "loess", size = .5) +
  coord_cartesian(ylim = c(0, 1)) +
  ylab("IPL-24 Integration Index") +
  xlab("Years of Residency") +
  ggtitle("") +
  theme_minimal() +
  theme(legend.position = "none",
        axis.title = element_text(size = 9),
        axis.text = element_text(size = 8)) 

mean(df$per240.5[df$yearsIn == 1])
mean(df$per240.5[df$yearsIn == 20])
mean(df$per240.5[df$yearsIn == 40])
mean(df$per240.5[df$yearsIn == 40]) - mean(df$per240.5[df$yearsIn == 20])
mean(df$per240.5[df$yearsIn == 20]) - mean(df$per240.5[df$yearsIn == 1])


ggsave(here::here("Draft/PNAS/draft/plots", "qrPlot12.pdf"), plot = QRplot12)
ggsave(here::here("Draft/PNAS/draft/plots", "qrPlot24.pdf"), plot = QRplot24)

